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Abstract 

A consistent treatment of both multiple scattering and small x quantum 
evolution effects on pair production in high energy pA collisions is feasi- 
ble in the framework of the Color Glass Condensate (CGC) [1]. We first 
discuss the properties of quark pair production in the classical effective 
theory where only multiple scattering effects are included. Explicit results 
are given for pair production as a function of the invariant mass of pairs, 
the pair momenta, the atomic mass number A and the quark mass. We 
relate the logarithms that appear in our formulation of pair production 
to logarithms that appear in the limit of coUinear factorization in QCD. 
Violations of k± factorization and medium modifications, as represented 
by the Cronin effect, are also investigated. We next consider how small 
X quantum evolution (shadowing) effects modify the results for pair pro- 
duction. In particular, we provide results for the rapidity distribution of 
pairs and the dependence of the Cronin effect on rapidity. We discuss the 
dependence of our results on the initial conditions for small x evolution 
and comment on its implications for pair production at RHIC and the 
LHC. 
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1 Introduction 



Proton-nucleus collisions provide a laboratory for studying the interaction of 
colored partons with an extended colored medium. Due to the combined effects 
of quantum nieehanieal coherence over extended longitudinal distances, and the 
QCD evolution of nuclear wave- functions, the parton from the projectile can 
probe both the density of color charges in the medium, and the multi-parton 
correlations that are intrinsic to high parton density components of the nuclear 
wave function. At high energies, the typical momentum transfer from partons 
in the medium to the probe is no longer soft and is characterized by a semi-hard 
scale^ 3> ^qcb- This scale, termed the saturation scale, is proportional to 
the density of partons in the transverse radius of the nucleus, and grows with 
energy [2, 3, 4]. Because the running of the coupling is controlled by this scale, 
asymptotic freedom tells us that the coupling of the colored partonic probe 
should be weak and will become weaker at higher energies. Therefore, with 
some effort, one can hope to learn about the properties of the medium - in the 
sense that one can compute reliably medium effects on the final state observables 
that are measured by experiments. 

The study of high parton density effects in QCD can be formulated, in weak 
coupling, as an effective field theory - the Color Glass Condensate (CGC) [5, 
6, 7, 8, 9, 10, 11, 12, 13]. The CGC has been widely applied to study gluon 
and quark production in proton-nucleus collisions- for a review, see ref. [14]. An 
attractive feature of the CGC effective theory is that one can quantify what 
one means by dilute or dense scatterers as a function of energy and mass num- 
ber [15]. What we mean by proton-nucleus collisions specifically, is a systematic 
expansion of amplitudes to lowest order in the ratio of the saturation momen- 
tum of the proton to the typical transverse momentum exchanged by the proton 
in the reaction {Qs,p/k_\_^p <C 1) and all orders in the ratio of the saturation mo- 
mentum of the nucleus relative to the momentum exchanged by the nucleus 
{Qs/k±.A)- At very high energies, or in very forward kinematics (in the frag- 
mentation region of the nucleus), the proton saturation scale can be large. In 
these kinematics, for fixed impact parameter, proton-nucleus collisions will be 
indistinguishable from nucleus- nucleus collisions. 

For gluon production in proton-nucleus collisions, the cross-section can be 
expressed in fcx -factor ized form as a product of unintegrated k±-dependent dis- 
tributions for the proton and the nucleus, convolved with the hard scattering 
matrix element [16, 17, 18, 19, 20, 15]. For a dilute projectile, such as a proton, 
the corresponding unintegrated gluon distribution at lowest order is a leading 
twist quantity which, integrated over kj_, gives the usiial leading log coUinear 
gluon distribution. In contrast, the unintegrated gluon distribution of the dense 
target, the nucleus, contains all twist contributions and has no analogue in the 
leading twist coUinear factorization formalism. In addition, this remarkable fac- 
torization is unlikely to hold beyond leading order in the expansion in powers 

^Throughout this paper, wc simply denote Qs the saturation niornentum of the nucleus. 
In the rare instances in which we need to introduce the proton saturation momentum, we will 
denote the latter Qs,p in order to distinguish from that of the nucleus. 
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of 0,,,p(xp)/fci,p [21, 22, 23]. 

It was shown in ref. [1] that /ej_ -factorization is broken expUcitly in pair 
production at leading order in proton-nucleus collisions. Novel 3-point and 4- 
point multi-parton ("all twist") correlation functions appear in the expression 
for the pair production cross-section. At large transverse momenta, k± ^ Qs, 
these expressions simplify [24] and the cross-section smoothly reduces to the 
fcj^-factorization result of Collins and Ellis [25] and Catani, Ciafaloni and Haut- 
mann [26]. The fc_L -factorization formalism of these authors has been used in 
several phenomenological studies of heavy quark production at collider ener- 
gies [27, 28, 29]. The magnitude of the breaking of fc_L-factorization for single 
quark production was quantified in ref. [30]. 

In this paper, we shall discuss in detail qualitative features of pair produc- 
tion that follow from the formalism developed in ref. [1]. A full treatment of 
multiple scattering and quantum evolution effects at high energies requires a 
computation of the previously mentioned 2-point, 3-point and 4-point correla- 
tion functions as a function of x, or of the rapidity Y (= ln(l/.T)). These can be 
determined in full generality (including all leading logarithms in x) by solving 
the Balitsky-JIMWLK equations for the small x evolution of multi-parton cor- 
relators [31, 32, 33, 34, 35, 36, 37, 38, 8, 9, 10, 11, 12, 13]. This would require 
an extensive numerical effort - only preliminary studies have been performed 
in this direction [39]. Nevertheless, a considerable deal can be learnt in certain 
limits. At large iV, and for large nuclei, the small x evolution of the 2-point cor- 
relators has a closed form expression called the Balitsky-Kovchegov (BK) equa- 
tion [37, 40, 41]. Numerical solutions of the Balitsky-Kovchegov equation are 
known in the fixed coupling [42, 43, 44, 45, 46, 47] and running coupling [39, 48] 
cases'^. In the large N and large A limit, the 3-point and 4-point correlators 
also simplify [54, 55, 56, 1] and their evolution in x can be expressed in terms of 
the solution of the BK equation. These arc computed numerically in this paper. 
We will also study heavy quark production in the McLerran-Venugopalan (MV) 
model [5, 6, 7], where high parton density effects contributing to multiple scat- 
tering are included but those arising from small x evolution arc not included. 
We are thus able to study quantitatively, in these MV and BK "mean field" 
limits, the interplay of multiple scattering and small x quantum evolution on 
pair production. 

The high parton density effects we discuss here include contributions cate- 
gorized as "higher twist" in the established language of coUinear factorization. 
These effects are difficult to treat in that framework and therefore a match- 
ing of the two formalisms is difficult. However, one can match the two for- 
malisms at large transverse momenta (specifically, when 4C }3^,M^) where 
fcj^-factorization formalism is recovered. Wc shall discuss in general how the 
logarithms that appear in the fcj^-factorization framework can be related to 
DGLAP [57, 58, 59, 60] logarithms that arise in coUinear factorization. Similar 
considerations arose in other specific processes [61, 62, 63, 64]. 

^ There are in addition several analytical studies which capture the key features of the BK 
equation [49, 50, 51, 52, 53]. 
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The formalism disciisscd here is vahd when at least one of the sources is 
dilute, as in pA collisions. When the dilute source becomes dense (either by 
increasing the energy of the proton or in a nucleus- nucleus collision), higher 
order corrections in Q,s,p/fc_L contribute significantly. Thus, in particular for 
nucleus-nucleus collisions, pair production has to be computed numerically; first 
results were obtained recently [65, 66]. 

We will not address phenomenological applications of our approach to pair 
production in pA collisions in this paper. There are interesting results from 
D-Au collisions at RHIC [67, 68] and in pA collisions at the SPS [69, 70, 71]. 
Some of this data has been studied previously in models based on the CGC 
approach [72, 73]. An interesting comparative study of several models of J/ if) 
production in heavy ion collisions and relevant references can be found in [74]. 
We will address the existing data and make predictions for future data on pair 
production at RHIC and LHC in a follow up to this work [75]. 

This paper is organized as follows. In section 2, we state the key results 
obtained in our previous derivation [1] of the pair production and single quark 
cross-sections. These were derived as a function of the momenta and rapidity 
of the quark and the anti-quark and we will re-express these in terms of the 
pair momenta, invariant mass and rapidity. In section 3, we shall discuss the 
properties of the multi-parton correlators and present results for their evolution 
as a function of energy in the large N and large A limit. In section 4, we 
will present results for the pair cross-sections. The results in this section are 
presented, in the CGC framework, for a) the McLerran-Venugopalan model and 
b) for the Balitsky-Kovchegov mean field model of the CGC. As mentioned 
previously, the former is a reasonable model in the kinematic region where 
multiple scattering effects are important and small x quantum evolution can 
be neglected. We discuss, in this model, pair distributions as a function of the 
pair momenta, rapidity and invariant mass. In particular, we show that the 
behavior of the spectra at large momenta can be understood analytically. The 
logarithms that appear in these expressions have the same origin as the collinear 
logs that appear in collinear factorization formalism of perturbative QCD. At 
smaller momenta, k±_ factorization is broken explicitly. We demonstrate the 
dependence of this breaking on the invariant mass and momenta of the pair. 
In the MV model, the saturation scale Qs depends on the size of the nucleus. 
We study the dependence of the pair cross-sections on Qs and on the quark 
mass. We next discuss the behavior of these quantities when small x quantum 
evolution is turned on. These shadowing effects are obtained by computing 
the rapidity dependence of correlators with the Balitsky-Kovchegov equation. 
We compute the rapidity distribution and the invariant mass and momentum 
distribution of pairs. In particular, we study the Cronin effect for quark pairs 
and its variation as a function of rapidity. We comment on the dependence 
of the small x quantum evolution on the large x initial conditions. Section 6 
summarizes our conclusions. 
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2 Quark cross-sections 



2.1 Generalities 

In the CGC formalism, the proton-nucleus collision is described as a collision 
of two classical fields originating from color sources representing the large x 
degrees of freedom in the proton and the nucleus. A schematic representation 
of the collision in this formalism is shown in fig. 1. As suggested by the figure, 



>A- 



Figure 1 : Cartoon representing gluon production in proton- nucleus collisions in 
the Color Glass Condensate framework. 

the color source distribution generating the classical field in each projectile has 
been evolved from some initial valence distribution at large x, in the leading 
logarithmic approximation in x, to the rapidity of interest in the collision. The 
gauge fields of gluons produced in the collision are determined by solving the 
Yang-Mills equations 

[D,,,F^''] = .r . (1) 

Here J'^ is the color current of the sources, which can expressed at leading order 
in the sources as 

J: = gS'^+Six-) pp^aixi_) + g5'^-5{x+)p,^a{xi_) , (2) 

where pp is the number density of "valence" partons in the proton moving in 
the -\-z direction at the speed of light. Likewise, p^ is the number density of 
"valence" partons in the nucleus moving in the opposite light cone direction. 
The previous two equations must be supplemented by a gauge fixing condition, 
and by the covariant conservation of the current : 

[D„,r]^Q. (3) 

The latter equation in general implies that eq. (2) for the current receives correc- 
tions that are of higher order in the sources pp and , because of the radiated 




5 



field. The solution of eqs. (1), (2) and (3) has been determined to all orders 
in both sources only numerically [22, 23, 76]. To lowest order in the proton 
source (as appropriate for a dilute proton source) and to all orders in the nu- 
clear source, analytical results are available and an explicit expression for the 
gauge field to this order, in Lorentz gauge, is given^ in ref. [15]. The amplitude 
for pair production to this order is obtained by computing the quark propa- 
gator in the background corresponding to this gauge field [1]. The diagrams 
corresponding to these insertions are shown in fig. 2. 




Figure 2: The two types of terms that enter in the QQ production amplitude 
in pA collisions. The gluon emitted by the proton can either emit the pair after 
the collision with the nucleus (left diagram) or before colliding with the nucleus 
(diagram on the right). The black blob denotes multiple scatterings resummed 
via Wilson lines. 

The probability for producing a single^ qq pair in the collision can be ex- 
pressed as 

''■i'-'''i-/(2|k/m:'^''''''''i'' 

where Adp{q,p) is the amputated time-ordered quark propagator in the pres- 
ence of the classical field. The argument [Pp,Pa\ indicates that this is the pro- 
duction probability in one particular configuration of the color sources. In order 
to turn this probability into a cross-section, one must first average over the initial 
classical sources pp and respectively with the weights ^ Wp[xp, pp], VF^ [x^ , p^] 
and subsequently integrate over all the impact parameters b, to obtain, 

a = j d?bj [Dpp] [DpJ Wp[xp, Pp] [x, , p,]Pi [pp, pj . (5) 

This formula incorporates both multiple scattering effects and quantum effects. 

The multiple scattering effects are in part included in the classical field (the 
solution of Yang-Mills equations), in part in the propagator of the quark in 

^This solution has been derived in [77] in the light-cone gauge of the proton, and in [61] in 
the gauge x~^A~ + A'^ = 0. 

^The probability of producing two pairs or more is parametrically of higher order in the 
proton source. Since we assume that pp is small, producing a single pair is thus the dominant 

process. 

''These weight functionals are normalized to ensure that their respective path integrals over 
the sources are unity. 
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this classical field, and also in the evolution with x of the source distributions. 
The quantum effects are included in the evolution of the weight functionals, 
Wp and , of the target and projectile with x. The arguments Xp and x^ 
denote the scale in x separating the largc-x static sources from the small-x dy- 
namical fields. In the McLerran-Venugopalan model, the functional that 
describes the distribution of color sources in the nucleus is a Gaussian in the 
color charge density ^ in [5, 6, 7]. Having a Gaussian distribution of sources, 
in our framework, is equivalent to the Glauber model of independent multi- 
ple scattering [15]. In general, this Gaussian is best interpreted as the initial 
condition for a non-trivial evolution of with x^. This evolution is 

described by a Wilson renormalization group equation - the JIMWLK equa- 
tion [31, 32, 33, 34, 35, 36, 37, 38, 8, 9, 10]. We will discuss evolution equations 
further in the following section. 



2.2 Differential pair cross-section 

The pair production cross-section, for quarks (anti-quarks) of momenta {p±) 
and rapidities yq = l/21n(g+/g~) {yp = 1/2 ln(p+/p~)), derived in [1] can be 
expressed as 



da 



(Pp^dPq^dypdyq 87t^{N^ - 1) 
trd 



"'l_L'*2_L 



{i+m)TqS-m)fT'\f ««"(fc2±; k^,k'^) 



trd 



(^+m)r,,-(2/-m)7°Tj7° + h.c. </>f «(fc2±; k^) 



+trd 



(6) 



where we denote ^ 



Tgg-(fei±,fe±) = 

Tg{ki±) 



2p+ [(g^ - + m2] + 2q+ [{q^ - fe± - fci± )2 + m2] 
_ ^Jp + g,fci^) 



(7) 



^This is true modulo a term proportional to the Cubie Casimir, that is parametrically 
suppressed relative to the Gaussian term at large A [78, 79]-sce eq. 20. 

*The momenta p and q of the produced particles have not been listed among the arguments 
of these objects to ensure the equations are more compa<;t. 
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Note that T'^^ = Tqq{ki±, k'j_) and is the well-known Lipatov vertex defined 
as 

C+{q, feii) = ^ + q+ ; C" {q, fei^) = ^ - q- ; {q, fei^) = -2k\ + q' , 

(8) 

with k2± = q± — ki±. 

<fp is the unintegrated gluon distribution of the proton and the various (j)^ 's 
arc unintegrated distributions describing the target. As stated previously, we 
are working to lowest order in the color charge density pp of the proton. The 
proton unintegrated distribution is given by ^ 



Mil-) ^ ^^-f^ / e*'"-- (p«(0)p«(x^)) . 
The corresponding nuclear "unintegrated distributions" are defined to be 

9 ^ Jx±,yj^ 



(9) 



xtv(Uix^)t''U^{y^)t''UbaiO) 

2_2r)2/2 r- 

(l)''^'i''(^lj_;k±,k'j_) — / gK-k±-y±+i±-{y±-y'±)+k'j^-(v'±-x'±)) 

9 N Jx'j^,yj_,y'j^ 

xiv(u{0)eu\y^)U{x'^)eu\y'S) ■ (10) 

The matrices U are path ordered exponentials along the light cone longitudi- 
nal extent corresponding to partonic configurations of the nucleus in the infinite 
momentum frame: 



U{x±) = V+ exp 



r + oo ^ 

/ dz+—^p^{z+,XA_)-T 

J -OO v_|_ 



(11) 



where the T° are the generators of the adjoint representation of SU (N) and 1^+ 
denotes a "time ordering" along the axis. The U have an identical definition, 
with the generators T" replaced by the generators t" in the fundamental repre- 
sentation of SU{N). We remind the reader that the expectation values ( • • • ) 



'in the MV model, the correlator of number densities is {pp{Q)pp{x±)) = p^(N^ 
1)5(2) (a;^ _ y^)^ where n'j^ = . Often, in the literature, the correlator of charge den- 

sities p,^ is used: /i^ = g'^l-i\- The unintegrated gluon distribution in eq. (9) is normalized 
such that the leading log gluon distribution in the proton satisfies 

r<3' 



1 

xGp{x,Q'^) = — J dl'j_(pp{l±) 
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here correspond to the averages over the sources pp in the case of ipp and in 
the case of the </>^ 's. 

In the definition of these distribution functions, the momentum l±_ is the total 
transverse momentum exchanged between the nucleus and the probe (either a 
gluon or the QQ state), and k±^k'^ are the momenta exchanged between the 
quark and the nucleus. These conventions can be visualized in the diagrams 




It is transparent that eq. (6) for the pair production cross-section is not in 
general fcj^-factorizable into a simple convolution of unintegrated parton distri- 
butions from the proton and the nucleus. While one can still factorize out the 
proton unintegrated distribution, the nucleus is now represented by the distribu- 
tions 0^'^, (/i'^'^, and 0^^''^'^, which are respectively 2-point, 3-point and 4-point 
correlators of Wilson lines in the color field of the nucleus. These multi-parton 
correlation functions contain all twists and are in general rather complicated. 
They however, in all generality, satisfy the sum rule, 

j cl)f'^Hlr,k^.k'^) - J ^Til±;k^) - ■ (13) 

k±,k'^ k± 

2.3 Leading twist limit 

It is illustrative to consider when one recovers fcj^-factorization. For Gaussian 
correlations, 

{p^,aix±_)p^,a'{x'^)) = 6aa'Sixj_ - x'^)^l\ , (14) 
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where /i^ = A/2'n:R\ cx A^f^. One obtains, in the leading twist approximation 
(where the path ordered exponentials are expanded to lowest order in p^/k'^), 



6i'««(fe2±;fei) 



TTi?' 



2 r 



''2X 



bf''^^k2±;k^,k'^)^7rR' 



^2± 



c 



X (27r)^ -f (5(fei - fe2±)5(feY - fe2±) + <5(fei)5(feY)) 



((5(fei - fe2±)5(fel) + (5(fci)5(fcl - fe2±)) 



(15) 



We assume here that the nuclei have a large uniform transverse area 7ri?^. The 
leading twist expressions for ^^'^^ and (pi^'i^ have a simple interpretation. At 
this order, the probe (gluon or qq pair) interacts with the nucleus by a single 
gluon exchange. The two delta functions in (ps^i^ correspond to the gluon being 
attached to the anti-quark line {S{k±) - there is no momentum flow from the 
nucleus to the quark line) or to the quark line {S{k2± — k±)- all the momentum 
from the nucleus flows on the quark line). An identical interpretation holds for 
the four terms of If one substitutes these leading twist approximations 

in cq. (6), one recovers the leading twist A;_L-factorized cross-section for pair 
production [25, 26, 24]. 



2.4 Kinematics 

We have to specify the multi-parton correlators ^^'^(f±), (t)''J'^{l±;k±), and 
(j)'ij''3^[l±;k±,k'j_) in order to compute the pair production cross-section. The 
following section will be devoted to a discussion of how one does this. Before 
proceeding, however, we will first discuss the kinematic variables in terms of 
which the pair production cross-section is specified. 

It is convenient to discuss the properties of the produced pairs in terms of the 
pair invariant mass , the pair rapidity Y and the pair transverse momentum 
P±. These are defined in terms of the momenta and rapidities of the quark and 
the anti-quark as 

P± = p± + q± 

2 \p-+q-) 2 \ujpe-yp +ojqe-yo ) 
= ujl+ujl + 2ujp ojg cosh {yp - Vg) - PI ■ (16) 

Here uig = \/q^+ and ojp = \/p^ + are the transverse masses of 
the quark and the anti-quark respectively, and yp = ln(p+/p~)/2 and yg = 
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\\i{q^ / q~) /2 arc their respective rapidities. We would like to express the cross- 
section in eq. (6) in terms of the pair kinematical parameters. Firstly, one notices 
that the cross-section in eq. (6) is expressed in terms of six variables while the 
pair invariants correspond to four variables. In other words, a given set of kine- 
matical parameters for the pair corresponds to a 2-dimensional manifold in the 
phase-space of the quark and anti-quark. One writes 



^ [ dq [ d(p — 

dM"^ d'^P± dY J J u}pU)q\ sinh (t/p - yg) \ d^px_dyp d'^q^^dyq 



(17) 

In the QQ cross-section that appears under the integral in eq. (17), the momenta 
q and p of the quark and anti-quark are given by 



= LM'i LAP.Ta Pom , (18) 
where q^^ and p'^^ are the momenta of the quark and anti-quark in the rest 

^10 



frame of the pair^ 



Qcm = I Y,9cos(/),gsm(/), Y — -m?-q^ 



Pcm= I Y.-9cos^;i,-g?sin(/),-y— -m2-g2 I . (19) 

To proceed from these momenta in the rest frame of the pair to the correspond- 
ing momenta in the laboratory frame in eq. (18), two Lorentz boosts must be 
applied, denoted by and L^. Assuming that the pair transverse momentum, 
P_L, is in the x direction^ ^, we first apply a boost in the x direction, with a ve- 
locity (3^ = \Pl_\/^JM'^ + P\_ - hence = + P\/M. At this point, the 
pair has a non-zero transverse momentum, but its (and hence its rapidity) is 
still zero. A second boost must be applied in the z direction, in order to bring 
the pair to the desired rapidity. The velocity of this boost in the z direction is 
f3z = tanh(y). The other factor in the integrand of eq. (17) is the Jacobian for 
the transformation pj^, q±,yp,yq P±,Y, M'^,q, (j). 

^"Wc arc here choosing to be positive and to be negative. However, the opposite 
choice is of course allowed as well. We don't need to consider it explicitly, and multiplying 
the final result by a factor 2 will be sufiicient. 

^^This choice is arbitrary, but has no influence on the result because the cross-section does 
not depend on the direction of fx ■ 
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3 Multi-parton correlators in the CGC 



3.1 Correlators in the MV model 



In the CGC approach, there is a procedure to compute the 2-, 3- and 4-point 

multi-parton correlators in cq. (10), in full generality, in the leading logarithmic 
approximation in x, starting from a given initial condition at some xq. However, 
no analytic solutions of these evolution equations are available and only prelim- 
inary work has been done in solving them numerically [39]. Fortunately, the 
multi-parton correlators simplify greatly in in the large N and large A asymp- 
totic limit of the theory. While this limit is truly asymptotic, we will assume 
that its domain of validity can be extended to finite A and N. 

In the kinematic domain \n{A^^^) <C <C A^^^ where = ln(l/a;^), the 
weight functional W^[x^, p^] has the form [5, 6, 7, 78], 



(-/ 



2ii 



^abc 



Pl{x)p\{x)p'^,{x) 



where (see eq. (14)) 



A 



27ri?2 



and K. 



A^N 



. (20) 



(21) 



The cubic term in the weight functional is parametrically suppressed by A^^^, 
and since we are working in the limit that A^^^ ^ 1, we will restrict ourselves 
in the rest of the discussion to the Gaussian term alone. 

The lower limit of the stated kinematical domain in rapidity follows from the 
requirement that a;^ <C which ensures that the small x partons in the 

nucleus couple coherently to all the color charges present along the z direction. 
The upper bound follows from the constraint that quantum corrections with a 
leading logarithm in x are small, namely, ln(l/a;^) <C 1. This condition com- 
bined with the large nucleus condition, a^A^^^ » 1, leads to the stated upper 
bound. In this kinematical domain in rapidity, for large nuclei, the model of 
Gaussian correlations of classical "valence" charges in the nuclear wave-function 
- the McLcrran-Vcnugopalan model - is valid. 

The MV model is thus a plausible model at moderately small values of x 
{x ~ 10~^), where small x quantum evolution effects are not yet large. It gives 
reasonable results for the initial conditions in heavy ion collisions at RHIC 
[22, 23, 80, 81, 82] and the Cronin effect in Deuteron-Gold collisions [62, 83, 84, 
85, 86, 87, 42], also at RHIC. In this Gaussian approximation, one can compute 
the 2-point, 3-point and 4-point correlators in closed form [38, 88, 89, 54, 55, 
56, 90, 1] - see Appendix A of ref. [1] for the explicit expressions for arbitrary 
N. We quote here the large N limits for these correlators because only these 



^^Note that the number density is related to the charge density by pA 
discussion of different conventions, see Ref. [22]. 



: gpA- For a 
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are used in the discussion of the energy evolution of correlators^^ : 



JV— >oo 



where we denote 



with 



Go(a;j_ - z±) = 



Go(y_L - z±) - Go{z± - u±) 



(22) 



(23) 



(24) 



tv(^Uix±)eu\y^)t'U^''\u^) 
tr (t''U'"'{x^)t'''u^'''''{u±) 



The other two correlators are obtained as special cases of the argument of the 
4-point function C 



C{x±,yj_;u±,u±) 
C{x±,x±;u±,u±) . (25) 



Thus in the large A'' limit, all the correlators that wc need in order to compute 
the cross-section in eq. (6) can be obtained from the quantity T{x± — yj_), 
defined by eqs. (23) and (24). 

When performing the Fourier transforms specified by eqs. (10) to obtain the 
(j)^ 's that enter in the cross-section, wc note that in the large N limit we have 



(26) 



This relation makes the sum-rule that relates the 3- and 4-point functions com- 
pletely obvious. Moreover, it makes the numerical integration simpler since it 
eliminates the integration over k' . For the purpose of numerical integrations, in 
the large N limit, the QQ cross-section can be efficiently written as 



da 



d^p^d^q^dypdyg Sn^N^ - 1) 



^{P± + ki± - k2±) 



k^ 



<pf'{k2±;k^){tTd (^+m)r,,-(2/-m)7°rt 7° 



+trd 



+trd 



{4+m)Tgg{i/-m)'y°Th° + h.c. 



-m)Tg(^-m)j'>Th° 



(27) 



i^The numerical evaluation of the QQ cross-section at finite N would be much more com- 
plicated and time consuming, because of the very complicated structure of the exact 4-point 
function (see [1] for the complete formula). 
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where wc have used the sum rule relating the 2- and 3-point functions to express 
the result in terms of the 3-point function alone. This expression is numerically 
very useful because exact cancellations that occur at small feij, among the three 
terms of eq. (6) are a consequence of the sum rules satisfied by the three (j)^ 's. 
This is difficult to ensure numerically and rewriting the cross-section by building 
in the sum rules, as in eq. (27), is an efficient way to achieve it. 

3.2 Energy evolution of correlators via the BK equation 

Recall that the previous results are limited to the regime where ag ln(l/x^) ^ 1. 
They therefore properly include multiple scattering effects which exist due to 
the large density of scatterers. However, they do not include the quantum effects 
arising from small- .x evolution, that generate leading twist shadowing. 

The assmnption of Gaussian correlations breaks down when asln(l/a;) ~ 
1 [91, 92]. The JIMWLK renormalization group equations [31, 32, 33, 34, 35, 
36, 37, 38, 8, 9, 10] incorporate, in the CGC framework, these large quantum 
corrections in the x evolution of VF[a;^, p^]. Consider for instance the correlator 
of two fundamental Wilson lines, {U{x±)W{yj_)), where the brackets denote 
the average with the nuclear weight functional W[x^, p^]. This correlator is 
directly proportional to the total qqA cross-section in deeply inelastic scatter- 
ing [88, 89] and appears as well in the single quark production cross-sections 
It satisfies the renormalization group equation, 

d\n{xo/x) \ I 2tt^ J {x± - z±y (z± - y^^y 

X (^Ntv[W{x^)U{y^)] - tr [U\x^)U{z^)] tr [W {z^)U{y^)]) . 

(28) 

This equation was first derived by Balitsky [37] and subsequently discussed in 
the CGC framework in Refs. [31, 32, 33, 34, 35, 36, 38, 8, 9, 10]. It cannot be 
solved in closed form since the 2-point function depends on the 4-point function,, 
and so on. However, in the large N and large A (a^A^^^ 3> 1) limit, the 
correlator of the product of traces of two pairs of Wilson lines factorizes into 
the product of the correlators of traces of pairs of Wilson lines: 

( tr [uHx^)U{z^)] tr [uHz^)Uiy^)] ) 

= {tv[U\x^)U{z^)]) {tT[U\z^)U{y^)]) , (29) 

to leading order in a 1/A^ expansion. With this factorization, eq. (28) becomes 
a closed form equation, called the Balitsky-Kovchegov (BK) equation [37, 41]. 

^■'Thc single quark cross-section can be obtained by integrating over the kinematic variables 
of one of the quarks in the pair [1, 30]. Using the sum rule in eq. (13), one can show that 
the correlators that contribute to single quark production are the 2-point correlators <p^^{l±) 
and (p^''{l±), and the 3-point correlator (p'^'^{l±;k±). The correlator <j>'^'> is the Fourier 
transform of {U{x±) U^iy±))- 
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While this equation has not been solved analytically (see however footnote 3), 
as discussed in the introduction, this equation has been solved numerically by 
several groups for both the fixed and running coupling cases [42, 43, 44, 45, 46, 
47, 39, 48]. For inclusive gluon production, where fcj^-factorization holds [20], 
the energy evolution of the two point correlator t/)^'^ has been studied using the 
BK evolution equation [42]. These studies show the same qualitative behavior as 
those seen in the remarkable measurements [93] in Deuteron-Gold collisions at 
RHIC of the disappearance of the Cronin peak with rapidity, and of the reversal 
of the centrality dependence of this effect from central to forward rapidities. 

These studies, with the BK equation, of the energy dependence of inclusive 
gluon production can be extended to quark pair production. This is because the 
factorization in eq. (29), essential to the derivation of the BK equation, requires 
that correlations between color charges be Gaussian. This may seem odd at first 
glance because we just argued that the MV model fails when as ln(l/a;) ~ 1. The 
resolution of this apparent paradox is that the correlations are not necessarily 
local anymore, namely, 

(p^,a(a;_L)p^,a'(y_L)) = ^ aa' IJ^\{x , X ^ - y . (30) 

The 2-, 3- and 4-point correlators are computed as in eqs. (22-25), with eq. (23) 
replaced by 

Ma '^C^' Vi- ~ '"-l) = / y (^z^Sr^ ftl {x, z± - r±) 

X (Go(y± - z±) - Go(m_l - z±)) 

X {Go{r^ - yj_) - Go{r^ - u^)) . (31) 

The a;-dependent l.h.s^ of this equation can be determined directly by solving 
the BK equation for {U{yj_) U^{u±)). Indeed, we have 

^f,lnx,y^-u^) = ~HU{y^)U\u^)) . (32) 

The r.h.s. of this equation is obtained by solving the BK equation in eqs. (28) 
and (29) with initial conditions given by the MV model. The l.h.s., thus deter- 
mined, specifies the values of the correlators in eq. (22) and eq. (25). Therefore, 
in the large A and large N limit, we can compute the energy dependence of all 
the correlators involved in pair production. This will enable us to study the 
effects of both multiple scattering and quantum evolution (shadowing) on the 
pair production cross-section. In the next section, we will discuss our results 
obtained in this approach. 

4 Results on pair production 

In this section, we will discuss results of our computation of the pair produc- 
tion cross-section in eq. (17) using the results for the correlators in the MV 
and BK models. The former includes multiple scattering effects but does not 
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include small x evolution. The latter includes both; it leads to a non-trivial 
A dependence, often called leading twist shadowing, that goes away only loga- 
rithmically with momentum. We would like to understand the dependence of 
pair production on the pair mass, the pair momentum, the quark mass and the 
rapidity of the pair. We would like to know how it depends on the saturation 
scale Qs- The saturation scale is a measure of the parton density in the system 
and depends on both x and A. In the first subsection, we will present results 
in the MV model to study the impact of the multiple scatterings alone on these 
quantities. In the second subsection, we will consider how small x quantum evo- 
lution a la BK modifies these results. As discussed previously, the MV results 
may be more relevant at central rapidities at RHIC; quantum evolution effects 
may be more relevant in forward Deuteron-Gold studies at RHIC and already 
for central rapidities at the LHC. 

4.1 MV model: multiple scattering effects 
4.1.1 M and dependence 




Figure 3: Pair cross-section (in arbitrary units) in the MV model, as a function 
of and M. 

In fig. 3, for illustration, we plot the cross-section for the quark pair produc- 
tion as a function of P_l and M . We have chosen here a quark mass of m=1.5 
GeV and the saturation scale to be Ql=2 GeV^. Unsurprisingly, it is peaked 

^^The dipole-hadron cross-section measured in DIS can be parameterized in terms of the 
saturation scale a la Golec-Biernat-Wusthoff [94, 95]. It is related to the scale /x^ in the 
MV model by the relation g'^fJ.\ = 4:1tQ1/Cf 1h{Qs/^qcd)' wh^^'s Cf is the Casimir in the 
fundamental representation. 
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at small values of P± and just above the threshold in the invariant mass M of 
the pairs. 

Examining the behavior of the differential pair cross-section at large M and 
fixed P±, or conversely at large P± and fixed M, we obtain the asymptotic 
forms, 

dN ln^(Af2) 

(33) 



P_[_ =const 

dN HPl) 

M=const ^ 

The l/P± and the l/Af^ behavior, in the stated limits is as one would expect in 
perturbative QCD (pQCD). Further, the logarithmic prefactors can be related 
to the collinear logarithms of pQCD which, at leading log order and beyond, 
are absorbed in the hadronic structure functions. 

The occurrence of these logarithms in the CGC framework and their rela- 
tion to the collinear logs of pQCD can be understood as follows. The CGC 
formulas obtained in [24] and [1] correspond to color sources from the proton 
and the nucleus that radiate gluons which fuse into a quark-antiquark pair. In 





Figure 4: Diagrams in the framework of collinear factorization obtained when 
one (top), or both (bottom), of the transverse momentum transfer to the pro- 
duced pair from the proton and nucleus becomes small. 

the limit where both the gluons from the proton and from the nucleus have a 
small transverse momentum, the CGC formulas reduce to the process gg — + QQ 
(at leading order) in the collinear factorization framework. Moreover, the inte- 
gration over the transverse momenta of these gluons produce logarithms that 
can be interpreted as the first power of the collinear logarithms resummed by 
the DGLAP equation. Because the large M limit is accessible at leading order 
in collinear factorization, it is natural to expect two powers of In(M^) in our 
formulas corresponding to the limit |fci_L|, |fc2_L| ~* 0. This is shown in the 
lower diagram of fig. 4. In contrast, the limit of large P± of the pair can only be 
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studied at next-to-leading order in the collinear factorization framework. In our 
formula, as shown in the upper diagram of fig. 4, this corresponds to the hmit 
where only one of the two gluons has a small transverse momentum (|fci_L| 0). 
Hence the single power of ln{P±). The relation of higher orders in the collinear 
factorization framework for pair production to the k± factorization framework 
was explored previously in refs. [26, 96]. 

4.1.2 Breaking of fcx-factorization 

Another important issue, previously studied for the production of single quarks 
in [30], is that of fcj^-factorization. This is illustrated in figure 5, where we 
display the cross-section versus P± for fixed M (left) and versus AI for fixed 
Pj_ (right). In these figures, we compare the exact values from eq. (17) (using 
eq. (6) and eq. (10) to the same expression in the fcj_-factorized approximation 
for pair production. The fc^-factorized result is obtained by replacing (j)^^'^^ and 




2 4 6 8 10 3 4 5 6 7 8 9 10 

Pi (GcV) M (GeV) 



Figure 5: Pair cross-section as a function of P± for fixed M (left). Pair cross- 
section as a function of M for fixed P± (right). The results in the fcj^-factorized 
approximation are shown with thin lines. 

(/>f in eq. (6) by 

(A^'99(fc2±;fei)fact = r/{k2±)(2nf [S{k^) + S{k^-k2±)] 



6f'99(fc2±;fci,fcl)fact = 0r(fc: 



'2A 



X (2^)4 {5{k^ - k2±)S{k'^ - k2±) + S{k^)5{k'^)) 

^ {S{k^-k2±)S{k'^)+S{k^)S{k'^-k2±)) 



2N 

(34) 

In other words, the limit of /c^-factorization is obtained from "nuclear distribu- 
tions" very similar to the leading twist formulas in eqs. (15), except that the 
leading twist 2-point function is replaced by the "all twist" 2-point function. 
This means that some of the rescattering corrections, but not all of them, can 
be included in the approximation of /c^-factorization. 



18 



The dependence of the ratio of the exact result to the fc_L-factorized result, 
as a function of P± and M, is nicely seen in the 3d plot of fig. 6. At any fixed 

ratio of pair Xsec, exact/fact, Q^, =2GeV", m=1.5GeV 




Figure 6: Ratio of the pair cross-section of the full result to the fc^-factorization 
result as a function of P± and M. 

M, for large P±, the exact cross-section and the fc^-factorized approximation 
become identical. (See also the left plot of figure 5.) This is because, in this 
limit, the quark and the anti-quark become coUinear with each of them having 
a very large transverse momentum. The quark-antiquark pair then scatters off 
the medium as a gluon would; as in the latter case, this leads to k± factorization. 

On the contrary, we observe in the right plot of figure 5 (fixed P± and 
large M), the exact cross-section and the fc_L-factorized approximation are not 
identical if the fixed value of the transverse momentum of the pair, P±, is of 
the order of Qs or smaller. This is because any pair configuration with a small 
total P± is very sensitive to rescatterings; even a small number of additional 
rescatterings, regardless of the pair invariant mass, may significantly change the 
transverse momentum of the pair. 

Turning now to fixed invariant mass and smaller transverse momentum 
P± < Qs, we note (see left panel of fig. 5) a qualitative change in the behav- 
ior of the cross-section due to multiple scattering, high parton density effects. 
In the fcj^-factorized approximation, the pair cross-section shows a bump at 
P± ^ Qs- Further it is suppressed relative to the exact result for P± < Qs- 
This suppression occurs because fc_L-factorization requires that only the quark 
or the anti-quark - and not both - scatters off the nucleus. The typical trans- 
verse momentum taken from the (dilute) proton is rather small; the transverse 
momentum of the pair is therefore approximately equal to the transverse mo- 
mentum exchanged during these scatterings on the nucleus-of order Qs- Thus, 
in the fcj^-factorized approximation, the pair is less likely to have a total mo- 
mentum P± smaller than Qs- In contrast, in the exact expression within the 
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MV model, both the quark and the anti-quark get multiply scattered in the tar- 
get nucleus. Therefore their net momentum can still be smaller than Qs- This 
explains the absence of a bump structure in the exact expression for the P± 
distribution of the quark-anti-quark pair. The size of the factorization breaking 
in this low P± region can be as large as a factor 3. 



violiition of fad for pairs, m=1.5GcV (l;u-gc N) violation of fad for pairs. m=4.5GeV (large N) 




4 6 8 10 12 14 16 18 20 8 10 12 14 16 18 20 22 24 

M (GcV) M (GcV) 



Figure 7: Ratio of the pair yields dN /dYdA'P from the exact formula and from 
the fcj_-factorized approximation as a function of M for to=1.5 GeV (left) and 
4.5 GeV (right). 

When we integrate the cross-section over the total momentum P ±^ , the mag- 
nitude of the factorization breaking, shown in fig. 7, is about 10% for m=1.5 
GeV and 2 GeV^. The latter is a typical value associated with central 

RHIC collisions. For the larger values of Qs that may be accessed in forward 
Deuteron-Gold collisions at RHIC and the LHC, the magnitude of the breaking 
is maximally 30%. Note that even for large invariant masses, the ratio returns to 
unity very slowly. This again is due to the fact that in the low P ±^ region, which 
contributes significantly to the integrated cross-section, the relative magnitude 
of violation of fc^-factorization is constant for large M . As for the case of single 
quark production [30], fig. 7 shows that the breaking of fc_L-factorization grows 
with increasing Qs while it is systematically weaker for larger quark masses. 

Finally, we observe in fig. 7 that, precisely at the pair threshold, the ratio 
of the exact to the fc^-factorized expression is unity. As is clear from the left 
panel of fig. 6, this ratio is not unity for any specific Pi but is a property of the 
distribution integrated over Pj_. This observation seems related to the following 
fact: at threshold, the quark and the antiquark are at rest in the rest frame 
of the pair. When boosted to the lab frame, they have the same momenta at 
threshold and are therefore coUinear. A pair made of a coUinear quark and 
antiquark is very similar to an octet gluon. The correlator of two Wilson lines 
satisfies a sum rule which ensures that the momentum integrated distributions 
are unchanged even if the momenta are redistributed by re-scattering. It is 
therefore plausible that the integrated pair distribution at threshold, in analogy 
to the gluon distribution, is insensitive to re-scatterings. 
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4.1.3 Nuclear size (Qs) dependence 

In the MV model, the nuclear size dependence arises only through the radius^^ 
R and the saturation scale Q^. In the MV model, the saturation scale 
is independent of the energy^^ and also depends on the atomic mass number 
roughly as A^^^. Alternately, the scale /Lt^ is used, which determines the number 
of color charges in the nucleus, per unit of transverse area. It is proportional to 
A^/^. The relation between the two scales is specified in footnote 15. 




Figure 8: Left: invariant mass spectrum of QQ pairs for different values of the 
saturation scale. Right: integrated quark yield as a function of /i^ ^ A^/"^, the 
scale specifying the number of color charges per unit transverse area. Points: 
numerical evaluation of the cross-section. Solid line: fit by a function in 

We first examine the integrated quark yield per unit rapidity (at fixed impact 
parameter) as a function of the nuclear size. This quantity is displayed in figure 
8, in differential form as a function of M on the left plot and integrated as a 
function of /x^ cx A^^^ on the right plot. Because one naively expects, for heavy 
quarks, a scaling with the number of binary collisions, namely A^^^, the curves 
in the left plot have been divided by a factor /i^ . The residual dependence on the 
saturation scale is therefore a departure from the binary scaling hypothesis. As 
one can see, when the saturation scale increases, the differential yield decreases, 
albeit in a very moderate fashion (note that the vertical axis for the left plot is a 
linear axis). When we integrate these functions over the invariant mass in order 
to obtain the total quark yield, we see that it is very close to a linear function 
of /x^ , confirming the fact that the departure from binary scaling is small. This 
departure is a small reduction of the yield compared to binary scaling, and 

^^This dependence on the nuclear size via the area of overlap between the two projectiles is 
trivial because translation invariance of all the multiparton correlators in the transverse plane 
is assumed. 

^^Albeit often, in saturation models, the Golec-Biernat-Wusthoff ansatz [Q^ = Qq (^)^, 
with A Si 0.3, and Qq = 1 GeV'^) is combined with the MV model. As discussed in section 3, 
this is consistent only if the Gaussian weight functional is non-local, as in the BK equation. 
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one can fit the dependence of dN/dY on //^ by a function that behaves hke 

Another quantity, whose dependence on the nuclear size is of phenomeno- 
logical interest, is the nuclear modification factor RpA- It is defined as the ratio 
between the yield in pA collisions to the yield in pp collisions, normalized by 
the number of binary collisions. 



R 



dN 

dM'^d'^P,dY 



■pA 



pA 



^1/3 dN 

dAPd'-T I dY 



(35) 



pp 



Indeed, this nuclear modification factor for pA collisions is essential for estab- 
lishing a baseline for "normal nuclear suppression" when looking for potential 
quark giuon plasma suppression effects on the yield of J/ip^s in nucleus-nucleus 
collisions. In the left plot of figure 9, we display this ratio as a function of P±, 



5 GcV, Q^- = 2 GeV" 




1.5 GcV. M = 3.1 GeV 



Px (GcV) 



Q,- (GcV-) 



Figure 9: Left: nuclear modification factor for the production of QQ pairs, as a 
function of their transverse momentum. Right: dependence of the modification 
factor at low P± as a function of Qs- The solid line is fit by a function that 
behaves as 1/Ql- 

for a given Qs and various fixed invariant masses. One can see here a behavior 
which is very similar to what was previously observed for gluon production in 
pA collisions (see [15] for instance): there is a suppression at low Pi and an 
enhancement at high P± , relative to perfect scaling with the number of binary 
collisions. 

The plot on the right of figure 9 shows how the value of the suppression 
at small momentum changes with varying saturation scale. We show this for 
an invariant mass just above the pair production threshold because this is the 
dominant kinematical domain contributing to the production of bound states 
(such as the J/'ip). We observe that the pair yield, for large Q^, behaves as 
l/Q^. Up to a logarithm (due to the small difference between and /Lt^), 
this corresponds to a scaling like 1/L where L is the longitudinal size of the 
nucleus. This 1/L dependence is well known from the super-penetration [97, 98] 
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of electron-positron pairs through metaUic foils (often called the "Chudakov 
effect") and was noted previously in the QCD case [99, 100, 55, 56], by looking 
at the propagation of QQ states through a random distribution of color fields. 
This result is in contrast to the Glauber like form exp{—C paL) that is often 
assumed, where a is the inelastic J/i/j-nucleon cross-section, p is the nuclear 
density, and C is constant number. Such an exponential form assumes successive 
independent collisions, and the decreasing exponential is thus the probability 
that the J /i}) has survived after propagating through a length L of nuclear 
matter. However, for the values of Qs probed at present energies it may be 
difficult to distinguish between the linear and exponential forms. 

1.2 
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Figure 10: Ratio of the pair cross-sections with the invariant mass M < 2Mo for 
"pp" and "pA." On the basis of the Color Evaporation Model, this is interpreted 
as the normal nuclear suppression factor of the charmonium. The parameters of 
the fitting curve are a = 1.51(6) GeV^ b = 2.65(16) GeV^ and a = 0.417(15). 

To observe how this pattern of suppression for producing QQ pairs translates 
into the nuclear J/^ suppression, we will consider the phenomenological "Color 
Evaporation" model (CEM)^'*. This model assumes that the yield of J/?/''s is 
obtained by integrating the yield of QQ pairs over the invariant mass of the 
pairs, from production threshold of a charm quark pair to the threshold for the 
production of a pair of D mesons, up to an overall constant prefactor. This can 
be expressed as 

dYd^Py_ CEM ^'"^ J^^, '^^ dYd^P^dA'P • 

If hadronization takes place outside of the nucleus (a legitimate assumption for 
J/ijj production at high energy), one can assume that the prefactor fj/^ is the 
same for pA and pp collisions. It therefore cancels out in the ratio RpA- In figure 

^*See [101] and references therein for a review on quarkonium production. 
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10, we represent the ratio RpA for J/tp production in the CEM as a function 
of the saturation momentum Qs ■ The sohd Hne represents a fit of this ratio by 
a power of Qs- One sees that the nuclear modification ratio is well reproduced 
by a power law, namely A""-^^ or 1/L° '*. We get a smaller power of 

Qs relative to figure 9 simply because we here integrate over Pi. The greater 
suppression at Pi = is attenuated by the enhancement present at higher Pi's. 
This pattern of Jf-ip suppression is not very sensitive to the precise mass of the 
c quark. 

4.1.4 Dependence on the quark mass 




Figure 11: Dots: pair yield as a function of the mass of the quarks. Solid line: 
fit by a function that behaves like ln(m)/TO'^. 

Finally, we shall discuss in the MV model, the dependence of the total QQ 
yield as a function of the mass of the quarks. This is plotted in figure 11 (per 
unit rapidity; as emphasized, there is no Y dependence in the MV model). 
One finds that, for quarks masses larger than the saturation scale, the yield 
behaves as ln(m)/m^, in agreement with what one expects from perturbation 
theory. The leading \/m? simply comes from integrating the l/PJ*; behavior 
from pi ~ m to oo. The additional factor of ln(m) arises from the intermediate 
Pi region, Qs < Pi < m where there is a deviation from the 1 /Pi behavior for 
masses below the saturation scale. The yield therefore fiattens out and does not 
grow as fast as 1/w? anymore. This implies that the limit of zero quark mass is 
less singular than it is at leading twist. In the latter case, the yield would keep 
growing as 1/m? all the way up into to the infrared region. 

^^The yield has been integrated over all J-'x before taking the ratio 
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4.2 Multiple scattering a la MV + quantum evolution a la 
BK 

4.2.1 Kinematics for small x evolution and initial conditions 

Our previous results, in the MV model, only included multiple scattering ef- 
fects on quark pair production in proton-nucleus collisions. Quantum evolution 
effects due to gluon bremsstrahlung are not included. At high energies, gluons 
with a very small longitudinal momentum fraction x may be probed in the pro- 
jectiles. It therefore becomes necessary to resum powers of asln{xo/x). These 
corrections cause leading twist shadowing of gluon distributions. Because the 
initial condition for small x quantum evolution is not the same in protons and 
nuclei, quantum evolution is yet another effect - in addition to the multiple 
scatterings already included in the MV model - that makes the cross-sections 
in pp and pA collisions differ. 

As discussed in section 3, we will study these quantum evolution effects 
in the framework of the Balitsky-Kovchegov (BK) equation. Wc noted there 
that the nuclear saturation scale, obtained from solutions of the BK equation, 
now depends on both A and the rapidity of the gluons probed in the nucleus. 
The former dependence comes from the initial condition, which is taken to be 
the McLerran-Venugopalan model, while the energy dependence follows from 
solutions of the BK equation. With the approximations for the correlators 
obtained by substituting the BK result for the l.h.s. of eq. (31) in eqs. (22-25), 
we are in a position to study how our previous multiple scattering results are 
modified by small x quantum evolution. 

In our kinematics, the -|- component of the gluon that comes from the proton 
is related to the -|- momentum of the pair by the relation 

k+=p+ + q+, (37) 

and similarly the — momentum taken from the nucleus is given by the relation 

K =p-+q- . (38) 

The longitudinal momentum fractions of the gluons probed in the proton and 
in the nucleus - respectively Xp and - are related to the -|- component of 
the momentum of the proton and to the — component of the momentum of a 
nucleon in the nucleus by : 



p+ ' p- 



(39) 



Using P+ = P~ = \/s/2, we can express these momentum fractions in terms of 
the transverse momentum of the quark and the antiquark and their rapidities 
(see the discussion following eq. (16)) as 

= . (40, 



25 



where Up = p\ + m? and uj"^ = q\ + . Alternatively, we can also express 
the momentum fractions in terms of the kincmatical parameters of the pair, 
Pl,M, Fas 

xp,, = y^^ e±^. (41) 

This latter form of Xp^A makes obvious the fact that when we integrate over q 
and (f) in eq. (17), the values of Xp and a;^ remain constant. This is very similar 
to gluon production, where the values of also fully determined by the 

final state^°. 

Our approach has an important caveat that we should discuss before pre- 
senting the results. We have argued elsewhere that the MV model, used here as 
the initial condition when solving the BK equation, is a good description of a 
nucleus at moderate value of x. This value is often estimated to be of order 
of xq = 10~^. By solving the BK equation, one obtains the value of the gluon 
distributions at any value of x smaller than this xq- However, in the calculation 
of the QQ cross-section, depending on the kinematics of the final state, we may 
need the unintegrated gluon distributions for values of Xp or x^ (or both) that 
are larger than xq. From eq. (41), we observe that this will happen in one of 
the projectiles if |y| becomes large, or in both of them at moderate Y and large 
Pt or M. Therefore, the MV initial condition at x = xq supplemented by BK 
evolution, is not sufficient to cover the entire kinematical range. One possibility 
would be to simply dismiss the points Pt,Y,M, where this problem occurs, on 
the grounds that these points are "contaminated" by large- .x physics. Ideally, 
one should try to constrain the necessary gluon distributions at large x from 
existing data. However, for the purposes of the qualitative discussion presented 
here, we follow an intermediate path: we make a simple ansatz for the gluon 
distributions that continues the MV initial condition at x > xq. 

Because the numerical evaluation of the pair cross-section is done with 
eq. (27), we need only to make such an ansatz for the proton unintegrated 
gluon distribution ipp and for the 3-point function (f>^^''^{k2j_', k±\x) in the nu- 
cleus. One generally requires, from quark counting rules, that for values of x 
close to unity, all gluon distributions vanish as {1 — x)'^. We adopt here a simple 
ansatz that implements this behavior; 

for x>xo, cf^f\k2r, k^\x) = (^f^) <^i'««"(fe2±; ki_\xo) , (42) 

where the denominator (1 — xq)^ ensures the continuity of </>^'''(fe2_L; k±\x) at 
x = xq- We make the same ansatz for (pp. Because the details of how these gluon 
distributions vanish when a; — > 1 are somewhat arbitrary, readers should note 

^'^The analogy between the two is not accidental and follows from the fact that all the 
produced particles are measured in both instances. In contrast, in single quark production, 
one integrates out the phase-space of the antiquark thereby allowing the values of Xp^^^ to 
vary. 

^"'This arbitrariness is similar to that of the initial condition Qq in DGLAP evolution. In 
both instances, asymptotic results are independent of the initial condition. 
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that results for the cross-section in kinematic regions dominated by large values 
of X (in one or both of the projectiles) contain a large systematic uncertainty. 

This uncertainty is reduced when the center of mass energy ^/s of the colli- 
sion increases because the values of Xp_^ decrease and the corresponding physics 
is less sensitive to the initial conditions. To ensure a significant region in the 
kinematical parameters P±,Y,M of the pair where the cross-section is not af- 
fected by this arbitrariness, all numerical calculations presented in the following 
subsections have been performed for the center of mass energy of the LHC, 
\/s = 8.5 TeV. At RHIC energies, if xq = 10~^, there is no region in phase- 
space which is not contaminated by the large x behavior of one (or both) of the 
projectiles. Therefore, a realistic description of RHIC results in this framework 
would require a less simplistic ansatz than eq. (42). Alternately, it may very 
well be the case that MV initial conditions may be applicable in large nuclei at 
larger values of xq. A systematic study is required but is beyond the scope of 
this paper. 

4.2.2 Rapidity dependence 

We shall now use the solution of the BK equation with the extrapolation done in 
eq. (42), to compute the rapidity dependence of the pair cross-section in eq. (17). 
In fig. 12, we show the rapidity distribution of cc pairs, for a fixed value of the 
invariant mass M = 3.1 GeV, and several values of the transverse momentum 
of the pair. The fact that the first derivative in x of the unintegratcd gluon 
distributions is not preserved by the ansatz of eq. (42) shows up in the cross- 
section as a discontinuity in the slope in Y of dN/dYdM^(pP±. On the plots, 
we have indicated the region in Y affected by this extrapolation by plotting 
only a solid line in those regions (their bomidaries can be easily determined 
from eq. (41)). The central region in rapidity - unaffected by this extrapolation 
- shrinks as one increases the transverse momentum of the pair, or its invariant 
mass. 

4.2.3 Cronin effect and shadowing 

To study deviations from the scaling with the number of binary collisions, one 
can compute the ratio RpA defined in eq. (35). In figure 13, we display this 
ratio as a function of Pi for different values of the rapidity (again the invariant 
mass is set to M = 3.1 GeV). On the left, the cross-section has been computed 
with the solution of the BK equation with fixed coupling, while on the right the 
BK equation was solved with a running coupling. There are small quantitative 
differences between fixed and running coupling, but the overall qualitative be- 
havior of RpA is unaffected. As one can see, the only region for which the ratio 
RpA goes over 1 is at the most negative rapidities. As the rapidity increases 

■^■^We have also tried an alternative large x ansatz of the form Ax" {I — x)**, where A and 
a are constrained to ensure that the unintegrated gluon distribution and its first derivative 
are continuous at x = xo- This ansatz in turn leads to odd behavior of the unintegrated 
distributions at large p± . 
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Figure 12: The rapidity distribution of quark-anti-quark pairs in pA collisions, 
for an invariant mass M = 3.1 GeV and various values of the transverse mo- 
mentum Pj^. The dots indicate the region in Y which is not sensitive to the 
extrapolation done at large x. (See eq. (42)). 
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Figure 13: Nuclear modification factor for QQ pairs at LHC center of mass 
energy, as a function of the transverse momentum of the pair. The pair invariant 
mass is just above the charm pair threshold {M = 3.1 GeV) and the rapidity is 
varied from —5 to -1-5. Left: the x dependence is obtained from the BK equation 
with fixed coupling [asNc/i^ — 0.2). Right: the x dependence is obtained from 
the BK equation with running coupling. Regions on the plot represented by 
only a solid line and no dots are sensitive to the large- a; extrapolation done in 
eq. (42). 
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above Y = 



—3, RpA decreases significantly, and remains below 1 at all Pi's. An 
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Figure 14: The rapidity dependence of the nuclear modification ration, at fixed 
M and for several values of Pi . The regions affected by the large x extrapolation 
are indicated by the absence of dots. 

alternative way to view the same information is to display RpA as a function of 
the rapidity Y, for fixed values of Pi, as in figure 14. (Note that for a limited 
region in Y, and for fixed Pi, it is proportional to the widely used variable xp-) 
On this plot, one sees a rapid decrease of RpA with Y. In our computation, 
this drop seems to start around Y — —3. Keep in mind however that the lower 
rapidity region is plagued by artifacts introduced by the large- a; extrapolation- 
the drop may in fact start earlier. Moreover, one sees on this plot that the 
decrease of RpA with Y is very rapid at negative Y, and slows down at positive 
rapidities. This is reminiscent of the very fast disappearance of the Cronin peak 
in the study of RpA for gluons via the BK equation [42, 20]. One also observes 
that the pattern of this suppression is almost independent of P± at small Pi, 
and that the Pi dependence becomes manifest only above Pi ~ 5 GeV. These 
results are qualitatively similar to those observed for the rapidity dependence 
of the nuclear modification factor at lower energies. 

5 Conclusions 

In this paper, we have explored the formalism developed in ref. [1] for pair 
production in proton nucleus collisions. We presented explicit results for two 
models, each of which corresponds to a particular limit of high parton densities 
realized in the Color Glass Condensate formalism of high energy QCD. In the 
first case, the McLerran-Venugopalan model, multiple scattering effects on pair 
production were studied in the absence of shadowing effects. We obtained re- 
sults in this model for the pair cross-sections as a function of the invariant mass, 
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the transverse momenta of the pairs, as weU as their dependence on nuclear size 
and the quark mass. For large transverse momenta or invariant masses, both the 
power law dependence on these scales, and the accompanying logarithms can 
be interpreted in terms of the leading kinematic behavior in the collinear fac- 
torization framework of perturbative QCD. We also studied, in the MV model, 
the extent of the violation of k± factorization. 

We next considered a model, where both multiple scattering and small x 
quantum evolution (shadowing) effects are included. In the limit of large Nc 
and large A, the energy evolution of the multi-parton correlators in pair cross- 
sections can be described in terms of the Balitsky-Kovchegov equation. We 
solved the BK equation numerically, for both fixed and running coupling, in 
order to determine the rapidity distribution of pairs as well as the evolution 
of the nuclear modification factor RpA as a function of rapidity. As in the 
gluon case, for fixed P±, one notes a rapid initial depletion of RpA with rapidity 
followed by a much slower depletion at larger rapidities. We emphasize that at 
RHIC energies our results may be sensitive to particular ansatze for unintegrated 
gluon distributions at large x; at LHC energies, this sensitivity is considerably 
weaker. 

The studies in this paper provide insight into the systematics of pair pro- 
duction and the dependence of the results on the quark masses, the system size 
and the center of mass energy. Detailed phenomenological studies of J/tjj and 
open charm production are underway and will be reported on in a foUowup to 
this work [75]. 

While this manuscript was being prepared for publication, ref. [102] ap- 
peared. The authors of ref. [102] correctly argue that the results of ref. [1] can 
be extended to include small x quantum evolution effects. This is precisely what 
is done here; quantum evolution effects are included by solving the BK equation 
for the unintegrated gluon distributions in eq. 10. Quantitative results from 
the solution of the BK equation for the evolution of pair cross-sections with 
energy/rapidity are presented in section 4.2. 
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